library(ggplot2)
## Warning: pacote 'ggplot2' foi compilado no R versão 4.4.3
library(dplyr)
## Warning: pacote 'dplyr' foi compilado no R versão 4.4.3
##
## Anexando pacote: 'dplyr'
## Os seguintes objetos são mascarados por 'package:stats':
##
## filter, lag
## Os seguintes objetos são mascarados por 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
## Warning: pacote 'tidyr' foi compilado no R versão 4.4.3
library(purrr)
library(tibble)
library(knitr)
library(bvartools)
## Warning: pacote 'bvartools' foi compilado no R versão 4.4.3
## Carregando pacotes exigidos: coda
## Carregando pacotes exigidos: Matrix
##
## Anexando pacote: 'Matrix'
## Os seguintes objetos são mascarados por 'package:tidyr':
##
## expand, pack, unpack
library(stringr)
set.seed(123)
# 1) Dados
data("us_macrodata") # contém Dp (inflação), u (desemprego), r (Fed Funds rates)
y <- cbind(Dp = us_macrodata[, "Dp"],
u = us_macrodata[, "u"],
r = us_macrodata[, "r"]) # ainda é 'ts' (mts) :contentReference[oaicite:1]{index=1}
# Estacionariza (exemplo) e define split com 'window()' para manter 'ts'
y_diff <- diff(y)
n_test <- 100
Ttot <- nrow(y_diff)
cut_id <- Ttot - n_test # última observação de TREINO
# datas correspondentes
t_all <- time(y_diff)
t_train_end <- t_all[cut_id]
t_test_start <- t_all[cut_id + 1]
y_train <- window(y_diff, end = t_train_end) # 'ts'
y_test <- window(y_diff, start = t_test_start) # 'ts'
# checagem rápida
stopifnot(inherits(y_train, "ts"), inherits(y_test, "ts"))
# índices na série estacionarizada
T_total_d <- NROW(y_diff)
stopifnot(n_test > 0, n_test < T_total_d)
idx_test <- (T_total_d - n_test + 1):T_total_d
idx_train <- 1:(T_total_d - n_test)
# Plot rápido das séries em nível (ou em taxa, no caso de Dp)
autoplot_y <- function(ts_obj, title) {
d <- as.data.frame(ts_obj)
d$time <- as.numeric(time(ts_obj))
d <- pivot_longer(d, -time, names_to="serie", values_to="valor")
ggplot(d, aes(time, valor)) +
geom_line() +
facet_wrap(~ serie, scales = "free_y", ncol = 1) +
labs(x = NULL, y = NULL, title = title)
}
autoplot_y(y, "Séries originais (Dp, u, r)")
Separar teste e treino:
# --- calcula ponto de corte (primeiro índice do teste) ---
ti <- time(y_diff)
split_idx <- min(idx_test)
split_time <- ti[split_idx]
# --- long data.frame para ggplot ---
df_long <- as.data.frame(y_diff) |>
mutate(time = ti) |>
pivot_longer(-time, names_to = "var", values_to = "value")
# --- plot com separação treino/teste ---
ggplot(df_long, aes(x = time, y = value)) +
# sombra do TESTE (vale para todos os painéis)
annotate("rect",
xmin = split_time, xmax = Inf, ymin = -Inf, ymax = Inf,
fill = "grey60", alpha = 0.12) +
# linha de corte
geom_vline(xintercept = split_time, linetype = "dashed", color = "grey40") +
# série
geom_line(color = "black", linewidth = 0.3) +
# facetas por variável com eixos livres
facet_wrap(~ var, ncol = 1, scales = "free_y") +
labs(
title = "Séries estacionarizadas (1ª diferença) — treino vs. teste",
subtitle = paste0("Treino: até ", sprintf("%.2f", split_time),
" | Teste: últimos ", length(idx_test), " períodos"),
x = NULL, y = NULL,
caption = "Faixa sombreada = período de teste"
) +
theme_minimal() +
theme(
strip.text = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
\[ \operatorname{vec}(Y)\sim\mathcal N\!\left((I_K\otimes X)\operatorname{vec}(\mathbf B),\;\Sigma\otimes I_T\right). \]
A log-verossimilhança \[ \ell(\mathbf B,\Sigma)\propto -\tfrac{T}{2}\log|\Sigma| -\tfrac{1}{2}\operatorname{tr}\!\big[\Sigma^{-1}(Y-X\mathbf B)'(Y-X\mathbf B)\big]. \]
Estimador OLS/MLE:
\[ \boxed{\;\widehat{\mathbf B}_{\text{MLE}}=(X'X)^{-1}X'Y.\;} \]
O MLE de \(\Sigma\) é
\[ \widehat{\Sigma}_{\text{MLE}}=\frac{(Y-X\widehat{\mathbf B})'(Y-X\widehat{\mathbf B})}{T}\quad \]
(ou com algum ajuste em T pra var.)
Para a \(k\)-ésima equação \(y_k = X\beta_k + \varepsilon_k\), com \(\varepsilon_k\sim\mathcal N(0,\sigma_k^2 I)\), um prior conjugado escalar:
\[ \beta_k\,|\,\sigma_k^2 \sim \mathcal N(\beta_{0k},\,\sigma_k^2 V_0),\qquad \sigma_k^2 \sim \mathcal{IG}\!\left(\tfrac{\nu_0}{2},\,\tfrac{s_0}{2}\right). \]
Então
\[ V_{n}=(V_0^{-1}+X'X)^{-1},\qquad \beta_{nk}=V_n\left(V_0^{-1}\beta_{0k}+X'y_k\right), \]
e, no limite a prior flat para \(V_0^{-1}\to 0\), temos a posterior \(\beta_{nk}\to(X'X)^{-1}X'y_k\) (OLS/MLE).
Como a prior afeta a posterior (para \(\mathbf B\)) A média posterior \(\mathbf B_n\) é uma média ponderada entre a informação amostral \(X'Y\) e a prior \(V_0^{-1}\mathbf B_0\), com precisão total \(V_0^{-1}+X'X\). Quando \(V_0\) é pequeno (grande \(V_0^{-1}\)), há maior shrinkage em direção a \(\mathbf B_0\) (para shrinkage em geral hyperprior=zero). Quando \(V_0\) é grande (pequeno \(V_0^{-1}\)), a prior pesa pouco e a posterior se aproxima do MLE.
\[ V_0^{-1}\to 0 \quad (\text{variância prior } V_0 \text{ muito grande}) \]
temos,
\[ V_n \to (X'X)^{-1},\qquad \mathbf B_n \to (X'X)^{-1}X'Y = \widehat{\mathbf B}_{\text{MLE}}. \]
Prior normal independente por coeficiente, com média típica \(m_{ij\ell}\) (geralmente \(m_{ii1}=1\) se a variável está em nível; caso contrário \(0\)) e variâncias que decaem com o lag e diferenciam “própria variável” vs. “cruzada”:
\[ b_{ij\ell}\;\sim\;\mathcal N\!\big(m_{ij\ell},\,\lambda_{ij\ell}\big),\qquad \lambda_{ij\ell} = \frac{\kappa_1}{\ell^{\kappa_3}}\times \begin{cases} \dfrac{\sigma_i^2}{\sigma_j^2}, & i=j,\\[6pt] \kappa_2\,\dfrac{\sigma_i^2}{\sigma_j^2}, & i\neq j~, \end{cases} \]
onde \(\kappa_1\) controla o encolhimento geral, \(\kappa_2\) o encolhimento de coeficientes cruzados e \(\kappa_3\) o decaimento com o lag \(\ell\).
A matriz de erros do VAR recebe uma prior conjugada do tipo Inverse–Wishart:
\[ \Sigma \;\sim\; \mathcal{IW}(S_0,\nu_0), \]
com hiperparâmetros \(S_0\) (matriz escala) e \(\nu_0\) (graus de liberdade). Em esquemas conjugados usuais, os coeficientes seguem uma Matrix–Normal condicionalmente a \(\Sigma\).
Prior de mistura spike-and-slab coeficiente a coeficiente: um indicador de inclusão escolhe entre variância “pequena” (spike) e “grande” (slab):
\[ \beta_k \,\big|\, \gamma_k \;\sim\; \mathcal N\!\big(0,\; \tau_{\gamma_k}^2\big), \qquad \gamma_k \sim \text{Bernoulli}(\pi),\qquad \tau_{\gamma_k}^2= \begin{cases} \tau_0^2,& \gamma_k=0 \ (\text{spike}),\\ \tau_1^2,& \gamma_k=1 \ (\text{slab}), \end{cases} \quad \text{com } \tau_0^2 \ll \tau_1^2. \]
set.seed(123)
# Especificações MCMC
p_lags <- 4
iters <- 2000
burnin <- 2000
# Base do modelo
base <- gen_var(y_train, p = p_lags, deterministic = "const",
iterations = iters, burnin = burnin)
# Priors Normal–Wishart com variâncias alvo: 1, 0.5, 0.25
# Atenção: v_i = 1/variância
K <- ncol(y_train)
priors <- list(
NW_var1 = add_priors(base,
coef = list(v_i = 1, v_i_det = 1),
sigma = list(df = K, scale = 0.0001)),
NW_var0.5 = add_priors(base,
coef = list(v_i = 2, v_i_det = 2),
sigma = list(df = K, scale = 0.0001)),
NW_var0.25 = add_priors(base,
coef = list(v_i = 4, v_i_det = 4),
sigma = list(df = K, scale = 0.0001))
)
# Prior SSVS (ajuste hiperparâmetros conforme preferência)
priors$SSVS03 <- add_priors(
base,
ssvs = list(inprior = 0.3, tau = c(0.1, 10), covar = FALSE, exclude_det = TRUE),
sigma = list(df = K, scale = 0.0001)
)
priors$SSVS05 <- add_priors(
base,
ssvs = list(inprior = 0.5, tau = c(0.1, 10), covar = FALSE, exclude_det = TRUE),
sigma = list(df = K, scale = 0.0001)
)
priors$SSVS08 <- add_priors(
base,
ssvs = list(inprior = 0.8, tau = c(0.1, 10), covar = FALSE, exclude_det = TRUE),
sigma = list(df = K, scale = 0.0001)
)
# Prior Minnesota (valores padrão razoáveis)
priors$Minn10 <- add_priors(
base,
coef = list(minnesota = list(kappa0 = 1, kappa1 = 1, kappa2 = NULL, kappa3 = 2)),
sigma = list(df = K, scale = 0.0001)
)
priors$Minn05 <- add_priors(
base,
coef = list(minnesota = list(kappa0 = 0.5, kappa1 = 0.5, kappa2 = NULL, kappa3 = 2)),
sigma = list(df = K, scale = 0.0001)
)
priors$Minn02 <- add_priors(
base,
coef = list(minnesota = list(kappa0 = 0.2, kappa1 = 0.2, kappa2 = NULL, kappa3 = 2)),
sigma = list(df = K, scale = 0.0001)
)
priors$Minn01 <- add_priors(
base,
coef = list(minnesota = list(kappa0 = 0.1, kappa1 = 0.1, kappa2 = NULL, kappa3 = 2)),
sigma = list(df = K, scale = 0.0001)
)
# Estima todos
mods <- lapply(priors, draw_posterior) # cada elemento vira um "bvar"
## Estimating model...
## Estimating model...
## Estimating model...
## Estimating model...
## Estimating model...
## Estimating model...
## Estimating model...
## Estimating model...
## Estimating model...
## Estimating model...
names(mods) <- names(priors)
# opcional: para usar summary(mods) como bvarlist
class(mods) <- c("bvarlist", "list")
names(mods)
## [1] "NW_var1" "NW_var0.5" "NW_var0.25" "SSVS03" "SSVS05"
## [6] "SSVS08" "Minn10" "Minn05" "Minn02" "Minn01"
Referências das funções usadas: gen_var,
add_priors, draw_posterior, Minnesota e
SSVS.
summary(bvar)retorna estatísticas das posteriores (pós burn-in). Vamos juntar tudo numa tabela e também fazer um histplot sobrepondo as densidades de um beta escolhido em cada modelo.
# ---- alvo: coeficientes de u defasada (L1..L4) na equação de u ----
eq_label <- "u" # equação alvo
var_label <- "u" # regressora alvo
max_lag <- 4 # quantos lags listar
# ---- helpers robustos ----
posterior_stats <- function(x) {
tibble::tibble(
Mean = mean(x),
SD = sd(x),
`2.5%` = stats::quantile(x, 0.025),
`50%` = stats::quantile(x, 0.50),
`97.5%`= stats::quantile(x, 0.975)
)
}
infer_dims <- function(obj){
A <- as.matrix(obj$A)
if (is.null(dim(A))) stop("Este 'bvar' não contém draws de A (A sem dimensão).")
S <- nrow(A); Pdim <- ncol(A)
# tenta inferir K por Sigma ou C
K <- NA_integer_
if (!is.null(obj$Sigma)) {
Sig <- as.matrix(obj$Sigma)
if (!is.null(dim(Sig))) {
k2 <- ncol(Sig)
Kc <- as.integer(round(sqrt(k2)))
if (Kc^2 == k2) K <- Kc
}
}
if (is.na(K) && !is.null(obj$C)) {
C <- as.matrix(obj$C)
if (!is.null(dim(C))) K <- ncol(C)
}
if (is.na(K) || K <= 0) stop("Não foi possível inferir K (Sigma/C ausentes).")
if ((Pdim %% (K*K)) != 0) stop(sprintf("ncol(A)=%d não é múltiplo de K^2=%d.", Pdim, K*K))
p <- as.integer(Pdim/(K*K))
vars <- colnames(obj$Y); if (is.null(vars)) vars <- colnames(obj$y)
if (is.null(vars)) vars <- paste0("y", seq_len(K))
list(K=K, p=p, vars=vars, S=S, A=A)
}
match_var <- function(label, vars){
idx <- match(label, vars)
if (is.na(idx)) idx <- match(tolower(label), tolower(vars))
idx
}
get_draws_idx <- function(obj, eq_idx, reg_var_idx, L){
d <- infer_dims(obj); K <- d$K; p <- d$p; A <- d$A; Kp <- K*p
if (L < 1 || L > p) stop(sprintf("Lag %d fora de 1..%d.", L, p))
col_idx <- (eq_idx - 1L)*Kp + (L - 1L)*K + reg_var_idx
as.numeric(A[, col_idx])
}
# ---- escolhe modelo de referência válido para obter p/vars ----
ref <- NULL
for (nm in names(mods)) {
tmp <- try(infer_dims(mods[[nm]]), silent = TRUE)
if (!inherits(tmp, "try-error")) { ref <- tmp; break }
}
if (is.null(ref)) stop("Nenhum modelo válido em 'mods' (sem draws de A).")
p_all <- ref$p
vars <- ref$vars
Ls <- seq_len(min(max_lag, p_all))
# ---- coleta estatísticas por modelo e por lag ----
stats_long <- purrr::imap_dfr(mods, function(m, nm){
dm <- try(infer_dims(m), silent = TRUE)
if (inherits(dm, "try-error")) return(NULL) # pula modelos problemáticos
eq_i <- match_var(eq_label, dm$vars)
reg_j <- match_var(var_label, dm$vars)
if (is.na(eq_i) || is.na(reg_j)) return(NULL)
dplyr::bind_rows(lapply(Ls, function(L){
draws <- get_draws_idx(m, eq_i, reg_j, L)
posterior_stats(draws) |>
dplyr::mutate(model = nm, lag = paste0("L", L), .before = 1)
}))
})
# versão longa (uma linha por modelo×lag)
stats_long <- stats_long |>
dplyr::arrange(lag, model)
# versão wide (cada lag vira um bloco de colunas)
stats_wide <- stats_long |>
dplyr::select(model, lag, Mean, SD, `2.5%`, `50%`, `97.5%`) |>
tidyr::pivot_wider(
names_from = lag,
values_from = c(Mean, SD, `2.5%`, `50%`, `97.5%`),
names_glue = "{.value}_{lag}"
) |>
dplyr::arrange(model) |>
dplyr::mutate(dplyr::across(-model, ~round(., 4)))
# se você tiver 'model_levels', força a ordem padrão dos modelos
if (exists("model_levels")) {
stats_wide <- stats_wide |>
dplyr::mutate(model = factor(model, levels = model_levels)) |>
dplyr::arrange(model) |>
dplyr::mutate(model = as.character(model))
}
stats_wide
## # A tibble: 10 × 21
## model Mean_L1 Mean_L2 Mean_L3 Mean_L4 SD_L1 SD_L2 SD_L3 SD_L4 `2.5%_L1`
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Minn01 -0.0296 0.0071 0.0053 -0.0523 0.0905 0.0109 0.0228 0.0716 -0.202
## 2 Minn02 -0.0199 0.0195 0.0133 -0.0593 0.108 0.0188 0.0408 0.0895 -0.227
## 3 Minn05 0.0053 0.0456 0.0148 -0.0814 0.125 0.0272 0.0749 0.111 -0.231
## 4 Minn10 0.0154 0.0581 0.0024 -0.0964 0.134 0.0311 0.0924 0.132 -0.243
## 5 NW_var… 0.0492 0.0607 0.0032 -0.0911 0.133 0.0333 0.104 0.135 -0.217
## 6 NW_var… 0.0465 0.0614 -0.0011 -0.102 0.131 0.033 0.103 0.140 -0.208
## 7 NW_var1 0.043 0.0633 -0.0094 -0.119 0.140 0.0329 0.104 0.148 -0.224
## 8 SSVS03 -0.0077 0.0429 0.0096 -0.0222 0.0798 0.0283 0.0671 0.0769 -0.166
## 9 SSVS05 -0.0073 0.0467 0.0074 -0.0235 0.0796 0.0298 0.0692 0.0767 -0.162
## 10 SSVS08 -0.0089 0.0543 0.0002 -0.0216 0.0809 0.0311 0.0715 0.079 -0.166
## # ℹ 11 more variables: `2.5%_L2` <dbl>, `2.5%_L3` <dbl>, `2.5%_L4` <dbl>,
## # `50%_L1` <dbl>, `50%_L2` <dbl>, `50%_L3` <dbl>, `50%_L4` <dbl>,
## # `97.5%_L1` <dbl>, `97.5%_L2` <dbl>, `97.5%_L3` <dbl>, `97.5%_L4` <dbl>
# --- níveis na ordem que deseja aparecer (usa a ordem da lista 'mods') ---
model_levels <- names(mods)
# --- ajuda pra extrair o número de cada família ---
num_or_na <- function(x) suppressWarnings(as.numeric(x))
# Separa modelos por família
nw <- model_levels[grepl("^NW_", model_levels)]
ssvs <- model_levels[grepl("^SSVS", model_levels, ignore.case = FALSE)]
minn <- model_levels[grepl("^Minn", model_levels, ignore.case = FALSE)]
# Extrai "variância" (ou proxy) de cada família para ordenar do MAIOR -> MENOR
# NW_varX: usamos o próprio X como variância alvo
nw_val <- num_or_na(sub("^NW_var", "", nw))
# SSVSxx: usamos a prob. de inclusão (ex.: 03 -> 0.3); menor = mais shrink
ssvs_val <- num_or_na(sub("^SSVS", "", ssvs)) / 100
# MinnXX: usamos kappa (ex.: 01 -> 0.1); menor = mais shrink
minn_val <- num_or_na(sub("^Minn", "", minn)) / 100
# Paletas (claro -> escuro); depois vamos REVERTER pra dar mais escuro ao menor valor
pal_blue <- colorRampPalette(c("#daebf9", "#9ecae1", "#3182bd", "#08306b"))
pal_green <- colorRampPalette(c("#d9f1d9", "#a1d99b", "#31a354", "#005a32"))
pal_red <- colorRampPalette(c("#fed6d9", "#fcae99", "#fb6a4a", "#cb181d", "#99000d"))
model_cols <- setNames(rep("#999999", length(model_levels)), model_levels) # default
# --- NW: menor variância -> mais ESCURO ---
if (length(nw)) {
ord_nw <- order(nw_val, decreasing = FALSE) # menor -> primeiro
cols_nw <- rev(pal_blue(length(nw))) # reverte: escuro para menores
model_cols[nw[ord_nw]] <- cols_nw
}
# --- SSVS: menor inprior -> mais shrink -> mais ESCURO ---
if (length(ssvs)) {
ord_ss <- order(ssvs_val, decreasing = FALSE)
cols_ss <- rev(pal_green(length(ssvs)))
model_cols[ssvs[ord_ss]] <- cols_ss
}
# --- Minnesota: menor kappa -> mais shrink -> mais ESCURO ---
if (length(minn)) {
ord_mn <- order(minn_val, decreasing = FALSE)
cols_mn <- rev(pal_red(length(minn)))
model_cols[minn[ord_mn]] <- cols_mn
}
# --- Use assim em QUALQUER ggplot em que você tenha
# >>> usa a paleta definida anteriormente
stopifnot(exists("model_cols"), exists("model_levels"))
# ---------- helpers ----------
infer_dims <- function(obj){
A <- as.matrix(obj$A); S <- nrow(A); Pdim <- ncol(A)
K <- NULL
if (!is.null(obj$Sigma)) {
k2 <- ncol(as.matrix(obj$Sigma)); Kc <- as.integer(round(sqrt(k2)))
if (Kc^2 == k2) K <- Kc
}
if (is.null(K) && !is.null(obj$C)) K <- ncol(as.matrix(obj$C))
if (is.null(K)) stop("Não foi possível inferir K.")
if (Pdim %% (K*K) != 0) stop("Dimensões de A incompatíveis com K.")
p <- as.integer(Pdim/(K*K))
vars <- colnames(obj$Y); if (is.null(vars)) vars <- colnames(obj$y)
if (is.null(vars)) vars <- paste0("y", seq_len(K))
list(K=K, p=p, vars=vars, S=S, A=A)
}
match_var <- function(label, vars){
idx <- match(label, vars)
if (is.na(idx)) idx <- match(tolower(label), tolower(vars))
idx
}
get_draws_idx <- function(obj, eq_idx, reg_var_idx, L){
d <- infer_dims(obj); K <- d$K; p <- d$p; A <- d$A; Kp <- K*p
if (L < 1 || L > p) stop(sprintf("Lag %d fora de 1..%d.", L, p))
col_idx <- (eq_idx - 1L)*Kp + (L - 1L)*K + reg_var_idx
as.numeric(A[, col_idx])
}
# ---------- constrói DF de densidades para uma equação ----------
build_dens_df <- function(eq_target, mods){
# referencial (pega o 1º modelo válido)
ref <- NULL
for (nm in names(mods)) {
try({ ref <- infer_dims(mods[[nm]]); break }, silent = TRUE)
}
if (is.null(ref)) stop("Nenhum modelo válido em 'mods'.")
K <- ref$K; p <- ref$p; vars <- ref$vars
eq_idx <- match_var(eq_target, vars)
if (is.na(eq_idx)) stop(sprintf("Equação '%s' não encontrada. Opções: %s",
eq_target, paste(vars, collapse=", ")))
rows <- list(); skipped <- character()
for (nm in names(mods)) {
res <- try({
dm <- infer_dims(mods[[nm]])
if (dm$K != K || dm$p != p) stop("Dimensões diferentes do modelo de referência.")
bind_rows(lapply(seq_len(K), function(vj){
bind_rows(lapply(seq_len(p), function(L){
tibble(
valor = get_draws_idx(mods[[nm]], eq_idx, vj, L),
modelo = nm,
reg_var = vars[vj],
L = L
)
}))
}))
}, silent = TRUE)
if (inherits(res, "try-error")) {
skipped <- c(skipped, nm)
} else {
rows[[nm]] <- res
}
}
out <- bind_rows(rows)
# >>> níveis/ordem iguais aos dos gráficos (usa model_levels)
out$modelo <- factor(out$modelo, levels = model_levels)
attr(out, "vars") <- vars
attr(out, "p") <- p
attr(out, "skipped") <- skipped
out
}
# ---------- plota 3x4 para UMA equação (linhas = regressoras, colunas = lags) ----------
plot_eq_grid <- function(eq_label){
dd <- build_dens_df(eq_label, mods)
vars <- attr(dd, "vars"); p <- attr(dd, "p")
dd <- dd %>%
mutate(
reg_var = factor(reg_var, levels = vars),
lag = factor(L, levels = 1:p, labels = paste0("L", 1:p))
)
ggplot(dd, aes(x = valor, color = modelo, fill = modelo)) +
geom_density(aes(y = after_stat(scaled)), alpha = 0.20, adjust = 1) +
facet_wrap(~ reg_var + lag, scales = "free", nrow = 3, ncol = 4) +
labs(
y = "densidade (normalizada: máx = 1)",
x = paste0(eq_label, " ~ regressor defasado"),
title = paste0("Posteriores dos coeficientes (BETAS) por equação: ", eq_label)
) +
# >>> aplica a paleta definida anteriormente
scale_color_manual(values = model_cols, limits = model_levels, drop = FALSE, name = NULL) +
scale_fill_manual(values = model_cols, limits = model_levels, drop = FALSE, name = NULL) +
theme_minimal() +
theme(legend.position = "bottom")
}
# ---------- gera os 3 gráficos (um por equação) ----------
eqs <- infer_dims(mods[[1]])$vars # ex.: c("Dp","u","r")
plots <- lapply(eqs, plot_eq_grid)
for (p in plots) print(p)
ariância preditiva marginal no horizonte \(h\) vem da diagonal de
\[ \Sigma_h \;=\; \sum_{i=0}^{h-1}\Psi_i\,\Sigma\,\Psi_i', \qquad \text{logo}\quad \sigma^2_{h,j}=\big[\Sigma_h\big]_{jj}. \]
\[ \boxed{\; \operatorname{RMSE}_j(h)\;=\;\sqrt{\frac{1}{O}\sum_{o=1}^{O} e_{o,h,\,j}^{\,2}} \;} \]
\[ \boxed{\; \operatorname{MAE}_j(h)\;=\;\frac{1}{O}\sum_{o=1}^{O} \big|e_{o,h,\,j}\big| \;} \]
Assumindo previsão Normal marginal para \(y_{o+h-1,\,j}\) com média \(\hat y_{o,h,\,j}\) e variância \(\sigma^2_{h,j}\),
\[ \boxed{\; \operatorname{LPD}_j(h)\;=\;\frac{1}{O}\sum_{o=1}^{O} \log\!\Big[\mathcal N\!\big(y_{o+h-1,\,j}\,\big|\,\hat y_{o,h,\,j},\,\sigma^2_{h,j}\big)\Big] \;=\;-\tfrac{1}{2}\log(2\pi\sigma^2_{h,j})-\frac{1}{2O}\sum_{o=1}^{O}\frac{e_{o,h,\,j}^{\,2}}{\sigma^2_{h,j}} \;} \]
tab <- tribble(
~`Símbolo`, ~Descrição,
"$h$", "Horizonte de previsão (passos à frente).",
"$H$", "Horizonte máximo avaliado.",
"$o$", "Índice da origem de previsão.",
"$O$", "Número total de origens de previsão consideradas.",
"$j$", "Índice da variável no vetor do VAR.",
"$K$", "Número de variáveis do VAR.",
"$t$", "Índice temporal na amostra.",
"$p$", "Ordem do VAR (número de defasagens).",
"$y_{o+h-1,\\,j}$", "Valor observado da variável $j$ no passo $h$ a partir da origem $o$.",
"$\\hat y_{o,h,\\,j}$", "Previsão (média preditiva) da variável $j$ no horizonte $h$ a partir da origem $o$.",
"$e_{o,h,\\,j}=y_{o+h-1,\\,j}-\\hat y_{o,h,\\,j}$", "Erro de previsão da variável $j$ no horizonte $h$ e origem $o$."
)
if (requireNamespace("kableExtra", quietly = TRUE)) {
kable(
tab,
format = ifelse(knitr::is_latex_output(), "latex", "html"),
booktabs = TRUE, escape = FALSE,
caption = "Notação e índices usados nas métricas."
) |>
kableExtra::kable_styling(full_width = FALSE, position = "left")
} else {
# fallback sem kableExtra
kable(
tab,
format = ifelse(knitr::is_latex_output(), "latex", "html"),
booktabs = knitr::is_latex_output(),
escape = FALSE,
caption = "Notação e índices usados nas métricas."
)
}
| Símbolo | Descrição |
|---|---|
| \(h\) | Horizonte de previsão (passos à frente). |
| \(H\) | Horizonte máximo avaliado. |
| \(o\) | Índice da origem de previsão. |
| \(O\) | Número total de origens de previsão consideradas. |
| \(j\) | Índice da variável no vetor do VAR. |
| \(K\) | Número de variáveis do VAR. |
| \(t\) | Índice temporal na amostra. |
| \(p\) | Ordem do VAR (número de defasagens). |
| \(y_{o+h-1,\,j}\) | Valor observado da variável \(j\) no passo \(h\) a partir da origem \(o\). |
| \(\hat y_{o,h,\,j}\) | Previsão (média preditiva) da variável \(j\) no horizonte \(h\) a partir da origem \(o\). |
| \(e_{o,h,\,j}=y_{o+h-1,\,j}-\hat y_{o,h,\,j}\) | Erro de previsão da variável \(j\) no horizonte \(h\) e origem \(o\). |
# --- helpers (usando seus objetos 'mods', 'y_diff', 'idx_test') ---
infer_params <- function(obj){
A <- as.matrix(obj$A); S <- nrow(A); Pdim <- ncol(A)
Sig <- as.matrix(obj$Sigma); K <- NULL
if (!is.null(Sig)) { k2 <- ncol(Sig); Kc <- as.integer(round(sqrt(k2))); if (Kc^2 == k2) K <- Kc }
if (is.null(K) && !is.null(obj$C)) K <- ncol(as.matrix(obj$C))
stopifnot(!is.null(K), Pdim %% (K*K) == 0)
p <- as.integer(Pdim/(K*K))
C <- if (!is.null(obj$C)) as.matrix(obj$C) else matrix(0, nrow=S, ncol=K)
list(K=K, p=p, S=S, A=A, C=C, Sigma=Sig)
}
posterior_means_ACSigma <- function(obj){
ip <- infer_params(obj); K <- ip$K; p <- ip$p; S <- ip$S
# A_bar
A_bar <- array(0, dim=c(K,K,p))
for (s in 1:S){
A_s <- matrix(ip$A[s,], nrow=K, byrow=TRUE)
for (L in 1:p) A_bar[,,L] <- A_bar[,,L] + A_s[, ((L-1)*K+1):(L*K), drop=FALSE]
}
A_bar <- A_bar / S
# C_bar
C_bar <- colMeans(ip$C)
# Sigma_bar
Sig_bar <- matrix(colMeans(ip$Sigma), nrow=K, byrow=TRUE)
list(A_bar=A_bar, C_bar=C_bar, Sigma_bar=Sig_bar, K=K, p=p)
}
# MA coeficientes Psi_i até H-1 (Ψ_0=I, Ψ_i = sum_{l=1..p} A_l Ψ_{i-l})
psi_list <- function(A_bar, H){
K <- dim(A_bar)[1]; p <- dim(A_bar)[3]
Psi <- vector("list", H); Psi[[1]] <- diag(K) # Ψ_0
for (h in 2:H){
acc <- matrix(0, K, K)
for (l in 1:min(p, h-1)) acc <- acc + A_bar[,,l] %*% Psi[[h-l]]
Psi[[h]] <- acc
}
Psi
}
# previsão dinâmica H-passos com A_bar/C_bar
forecast_path_mean <- function(A_bar, C_bar, Y, t0, H){
K <- ncol(Y); p <- dim(A_bar)[3]
L <- matrix(NA_real_, nrow=p, ncol=K); for (ell in 1:p) L[ell,] <- Y[t0-ell,]
Yhat <- matrix(NA_real_, nrow=H, ncol=K)
for (h in 1:H){
y_new <- C_bar
for (LAG in 1:p) y_new <- y_new + drop(A_bar[,,LAG] %*% L[LAG,])
y_new <- as.numeric(y_new)
Yhat[h,] <- y_new
if (p>1) L <- rbind(matrix(y_new,1), L[1:(p-1),,drop=FALSE]) else L[1,] <- y_new
}
Yhat
}
metrics_by_h <- function(obj, y_all, idx_test, H=8){
pm <- posterior_means_ACSigma(obj)
A_bar <- pm$A_bar; C_bar <- pm$C_bar; Sig_bar <- pm$Sigma_bar
K <- pm$K; p <- pm$p; Y <- as.matrix(y_all); Tn <- nrow(Y)
# Índices válidos
idx <- if (is.ts(y_all)) {
pos <- seq_len(Tn); names(pos) <- as.character(time(y_all))
ii <- if (is.numeric(idx_test)) idx_test else match(idx_test, names(pos)); as.integer(ii)
} else as.integer(idx_test)
idx <- idx[!is.na(idx) & idx>p & (idx+H-1)<=Tn]
stopifnot(length(idx)>0)
# Prepara Ψ e Σ_h = sum_{i=0..h-1} Ψ_i Σ Ψ_i'
Psis <- psi_list(A_bar, H)
Sigma_h <- lapply(1:H, function(h){
Sacc <- matrix(0, K, K)
for (i in 1:h) Sacc <- Sacc + Psis[[i]] %*% Sig_bar %*% t(Psis[[i]])
Sacc
})
err2 <- lapply(1:H, function(...) matrix(0.0, nrow=length(idx), ncol=K))
absE <- lapply(1:H, function(...) matrix(0.0, nrow=length(idx), ncol=K))
lpdf <- lapply(1:H, function(...) matrix(0.0, nrow=length(idx), ncol=K))
for (o in seq_along(idx)){
t0 <- idx[o]
Yhat <- forecast_path_mean(A_bar, C_bar, Y, t0, H) # H x K
for (h in 1:H){
y_true <- Y[t0+h-1,]
e <- y_true - Yhat[h,]
err2[[h]][o,] <- e^2
absE[[h]][o,] <- abs(e)
sdv <- sqrt(pmax(diag(Sigma_h[[h]]), .Machine$double.eps))
lpdf[[h]][o,] <- dnorm(y_true, mean=Yhat[h,], sd=sdv, log=TRUE)
}
}
cn <- colnames(Y); if (is.null(cn)) cn <- paste0("y", seq_len(K))
bind_rows(lapply(1:H, function(h){
tibble(
horizon = h,
var = cn,
RMSE = sqrt(colMeans(err2[[h]])),
MAE = colMeans(absE[[h]]),
LPD = colMeans(lpdf[[h]]) # LPD marginal (por variável) médio nas origens
)
}))
}
# ---- roda para todos os modelos, H=8 ----
H <- 4
metrics_h <- purrr::imap_dfr(
mods, ~ metrics_by_h(.x, y_diff, idx_test, H) %>% mutate(model = .y)
) %>% dplyr::relocate(model, var, horizon)
metrics_h
## # A tibble: 120 × 6
## model var horizon RMSE MAE LPD
## <chr> <chr> <int> <dbl> <dbl> <dbl>
## 1 NW_var1 Dp 1 1.30 0.869 -7.24
## 2 NW_var1 u 1 0.385 0.263 -0.582
## 3 NW_var1 r 1 0.684 0.532 -1.22
## 4 NW_var1 Dp 2 0.946 0.629 -1.42
## 5 NW_var1 u 2 0.356 0.255 -0.386
## 6 NW_var1 r 2 0.627 0.482 -1.19
## 7 NW_var1 Dp 3 1.30 0.907 -1.68
## 8 NW_var1 u 3 0.505 0.373 -0.821
## 9 NW_var1 r 3 0.602 0.448 -1.19
## 10 NW_var1 Dp 4 0.833 0.611 -1.53
## # ℹ 110 more rows
# usa a paleta definida anteriormente
stopifnot(exists("model_cols"), exists("model_levels"))
stopifnot(all(c("model","var","horizon","RMSE","MAE","LPD") %in% names(metrics_h)))
var_levels <- unique(metrics_h$var)
Hmax <- max(metrics_h$horizon, na.rm = TRUE)
metrics_h_long <- metrics_h %>%
pivot_longer(c(RMSE, MAE, LPD), names_to = "metric", values_to = "value") %>%
mutate(
metric = factor(metric, levels = c("RMSE","MAE","LPD")),
var = factor(var, levels = var_levels),
# >>> força os níveis/ordem dos modelos para bater com a paleta
model = factor(model, levels = model_levels),
horizon_f = factor(horizon, levels = seq_len(Hmax), labels = paste0("h", seq_len(Hmax))),
facet_id = paste0(as.character(metric), " · ", as.character(var)),
label = ifelse(metric == "LPD", sprintf("%.2f", value), sprintf("%.3f", value)),
vjust_lab = ifelse(value < 0, 1.2, -0.25)
) %>%
group_by(facet_id) %>% # mesma escala por LINHA (todas as colunas)
mutate(
y_min_row = min(value, na.rm = TRUE) - 0.05,
y_max_row = max(value, na.rm = TRUE) + 0.05
) %>%
ungroup()
# ---- controles visuais ----
bar_lwd <- 4 # espessura da “barra” (segmento)
gap <- 1.30 # >1 aumenta o espaço entre modelos
# posição X por modelo (constante em todas as colunas)
M <- nlevels(metrics_h_long$model)
metrics_h_long <- metrics_h_long %>%
mutate(x_m = as.integer(model) * gap)
x_breaks <- gap * seq_len(M)
x_labels <- levels(metrics_h_long$model)
# ordem das 9 linhas: 3 RMSE, 3 MAE, 3 LPD
facet_levels <- c(paste0("RMSE · ", var_levels),
paste0("MAE · " , var_levels),
paste0("LPD · " , var_levels))
metrics_h_long$facet_id <- factor(metrics_h_long$facet_id, levels = facet_levels)
ggplot(metrics_h_long, aes(x = x_m, y = value, color = model)) +
geom_segment(aes(xend = x_m, y = y_min_row, yend = value),
linewidth = bar_lwd, lineend = "butt") +
geom_point(size = 2) +
geom_text(aes(label = label, vjust = vjust_lab), color = "black", size = 4) +
geom_blank(aes(y = y_min_row)) + # fixa limites da linha
geom_blank(aes(y = y_max_row)) +
facet_grid(facet_id ~ horizon_f, scales = "free_y") + # 9 linhas x 8 colunas
scale_x_continuous(
breaks = x_breaks,
labels = x_labels,
expand = expansion(add = c(gap*0.6, gap*0.6))
) +
# >>> aplica a paleta customizada (cores por modelo)
scale_color_manual(values = model_cols, limits = model_levels, drop = FALSE, name = NULL) +
scale_fill_manual(values = model_cols, limits = model_levels, drop = FALSE, name = NULL) +
labs(
x = "Modelos (em cada horizonte)", y = NULL,
title = "RMSE, MAE e LPD — 9 linhas × 8 colunas",
subtitle = "Colunas = horizontes; mesma escala Y por linha; paleta NW=azuis, SSVS=verdes, Minn=vermelhos"
) +
theme_minimal() +
theme(
legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
strip.text.y.left = element_text(angle = 0)
)
# precisa da paleta definida antes
stopifnot(exists("model_cols"), exists("model_levels"))
stopifnot(all(c("model","var","horizon","RMSE","MAE","LPD") %in% names(metrics_h)))
# helpers p/ ordenar dentro de cada facet
reorder_within <- function(x, by, within, fun = mean, sep = "___") {
x_within <- paste(x, within, sep = sep)
stats::reorder(x_within, by, FUN = fun)
}
scale_x_reordered <- function(..., sep = "___") {
scale_x_discrete(labels = function(x) gsub(paste0(sep, ".*$"), "", x), ...)
}
# média sobre variáveis e horizontes
metrics_mean <- metrics_h %>%
group_by(model) %>%
summarise(
RMSE = mean(RMSE, na.rm = TRUE),
MAE = mean(MAE, na.rm = TRUE),
LPD = mean(LPD, na.rm = TRUE),
.groups = "drop"
) %>%
pivot_longer(c(RMSE, MAE, LPD), names_to = "metric", values_to = "value") %>%
mutate(
metric = factor(metric, levels = c("RMSE","MAE","LPD")),
# força níveis dos modelos para casar com a paleta
model = factor(model, levels = model_levels),
label = ifelse(metric == "LPD", sprintf("%.2f", value), sprintf("%.3f", value)),
vjust_lab = ifelse(value < 0, 1.2, -0.25),
model_ord = reorder_within(model, value, metric)
) %>%
group_by(metric) %>%
mutate(
y0 = min(value, na.rm = TRUE) - 0.05,
mid = (value + y0) / 2,
height = pmax(abs(value - y0), .Machine$double.eps)
) %>%
ungroup()
ggplot(metrics_mean, aes(x = model_ord, fill = model)) +
geom_tile(aes(y = mid, height = height), width = 0.7) +
geom_text(aes(y = value, label = label, vjust = vjust_lab), color = "black", size = 3) +
facet_wrap(~ metric, scales = "free", nrow = 1) +
scale_x_reordered() +
# aplica a paleta
scale_fill_manual(values = model_cols, limits = model_levels, drop = FALSE, name = NULL) +
labs(
x = NULL, y = NULL,
title = "Média das métricas por modelo",
subtitle = "RMSE, MAE e LPD (cores conforme famílias/prior)"
) +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.05))) +
theme_minimal() +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)
)
par(mfrow = c(4, 4), mar = c(3,3,2,1))
nh <- 12 # horizonte
plot_irf_one <- function(obj, title) {
ir <- irf(obj, impulse = "u", response = "Dp", n.ahead = nh, type = "gir")
plot(ir, main = title)
}
model_names <- names(mods)
for (nm in model_names) {
plot_irf_one(mods[[nm]], nm)
}
# Se sobrar um painel, deixe em branco